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Abstract 

Lienard equations, x + ef(x)x + x = 0, with f{x) an even continuous function are consid- 
ered. In the weakly nonlinear regime (e — ► 0), the number and a (D(e ) approximation of 
the amplitude of limit cycles present in this type of systems, can be obtained by applying 
a methodology recently proposed by the authors [Lopez-Ruiz R, Lopez JL. Bifurcation 
curves of limit cycles in some Lienard systems. Int J Bifurcat Chaos 2000; 10:971-980]. In 
the present work, that method is carried forward to higher orders in e and is embedded 
in a general recursive algorithm capable to approximate the form of the limit cycles and 
to correct their amplitudes as an expansion in powers of e. Several examples showing the 
application of this scheme are given. 
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1 Introduction 



Self-sustained periodic oscillations are common in nature pQ. They are also important in 
engineering applications [2] • The calculation of the number and amplitude of such different 
periodic motions (limit cycles) taking place in an oscillating system is an unsolved problem. 
This question constitutes the second part of Hilbert's Sixteenth Problem jHJ Q] when we 
are restricted to two-dimensional autonomous systems of the form: 

x = P n {x,y), 

y = Qn(x,y), (1) 

where x(t) = dx(t)/dt, y(t) = dy(t)/dt and (P n , Q n ) are polynomials of degree n with real 
coefficients. Although it has been proved that the number of limit cycles in systems of type 
(fl|) is finite 0IB], the determination of the maximal number H(n) of such solutions for a 
given degree n is far away of being known. Even for n = 2, H (2) is still not determined 
It has been verified with different examples, for instance, that H{2) > 4, H(3) > 11, 
H(i) > 15 and H(5) > 23 but the exact values of H(n) for these cases El El EH] are 
unknown. 

The van der Pol oscillator x + e(x 2 — l)x + x = 0, is an example of system (Q) that has 
been well studied. In this case, Ps(x,y) = y and Q3(x,y) = —e{x 2 — l)y — x. It displays 
a limit cycle whose uniqueness and non-algebraicity has been shown for the whole range 
of the parameter e |llj . Its behaviour runs from near-harmonic oscillations for e close to 
zero (e — ► 0) to relaxation oscillations when e tends to infinity (e — > oo), making it a good 
model for many practical situations |12j . 

A generalization of the van der Pol oscillator is the Lienard equation, 

x + ef(x)x + x = 0, (2) 

with e a real parameter and f(x) any real function. When f(x) is a polynomial of degree 
= 2n + 1 or 2n this equation takes the form (1) with Pn+i(x, y) = y and Qn-\-i{x, y) = 
—ef(x)y — x. It has been conjectured by Lins, Melo and Pugh (LMP-conjecture) that the 
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maximum number of limit cycles allowed is just n ^SI- It is true if N = 2, or N = 3 or if 
f(x) is even and N = 4 |13| I14j. Also, there are strong arguments for claiming its truth 
in the strongly nonlinear regime (e — > oo) when f{x) is an even polynomial There 
are no general results about the limit cycles when fix) is a polynomial of degree greater 
than 5 neither, in general, when f(x) is an arbitrary real function 

The calculation of the form of limit cycles is another difficult task associated to this 
problem. As far as we know, there is no a general methodology documented in the litera- 
ture on this subject. The different schemes used to calculate the number and amplitude of 
the periodic motions in the weak nonlinear regime, namely the averaging and perturbation 
methods 117j , could also be applied to find the approximate wave form of those orbits in 
the time or frequency domain ^Hj. In general, the insistence in performing these calcula- 
tions in phase space with the time variable being explicit can complicate the achievement 
of the objective. 

In this work, we undertake the task of calculating the form of limit cycles in phase 
space for the weakly nonlinear regime of Lienard equations. The method exploited in 
|19| I2()j to find the number and a first order approximation to the amplitude of limit 
cycles when | e |<C 1 is recalled in Section 2. The embedding of this method in a general 
recursive framework for calculating the form and amplitude of those periodic solutions is 
presented in Section 3. The general algorithm is rewritten in Section 4. Some illustrative 
examples are given in Section 5. Finally, we present our conclusions. 



2 Integral equation for the limit cycles 

In order to study the limit cycles of equation (J2J) with the time variable being implicit, it 
is convenient to rewrite it in the coordinates (x, x) = (x, y) in the plane, with x(t) = y(x) 
and x(t) = y(x)y'(x) (where y'{x) = dy/dx): 

yy + ef(x)y + x = 0. (3) 
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A limit cycle C\ = (x,y±(x)) of equation © has a positive branch y+(x) > and a 
negative branch y~(x) < 0. They cut the x-axis in two points (— a_,0) and (a+,0) with 
a_,a+ > because the origin (0,0) is the only fixed point of Eq. (J3J). Then every limit 
cycle C\ solution of Eq. (j2J) encloses the origin and the oscillation x runs in the interval 
— a_ < x < a_|_. The amplitudes of oscillation a_,a + identify the limit cycle. The result is 
a nested set of closed curves that defines the qualitative distribution of the integral curves 
in the plane (x,y). The stability of the limit cycles is alternate. For a given stable limit 
cycle, the two neighbouring limit cycles, the closest one in its interior and the closest one 
in its exterior, are unstable, and viceversa. 

When f(x) is an even function, the inversion symmetry (x,y) <-> (— x, —y) of Eq. © 
implies y+(x) = — y~(— x) and then a\ = 0,2 = a. Therefore we can restrict ourselves to the 
positive branches of the limit cycles with —a < x < a. In this case, the ampli- 

tude a identifies the limit cycle. The parameter inversion symmetry (e, x, y) «-> (— e, x, —y) 
implies that if C\ = (x,y±(x)) is a limit cycle for a given e, then C\ = (x, —y^(x)) is a 
limit cycle for — e. In consequence, the amplitude a of the limit cycles in these Lienard 
systems is an even function of e. Moreover if C/ is stable (or unstable) then Ci is unstable 
(or stable, respectively). Therefore it is enough to consider the limit cycles when e > for 
obtaining all the periodic solutions. (The limit cycles for a given e < are obtained from 
a reflection over the x-axis of those limit cycles obtained for e > 0). 

Another global property of a limit cycle of Eq. can be derived from the fact that 
the mechanical energy E = (x 2 + y 2 )/2 is conserved in a half oscillation: 



By Eq. (j3J), we have 4^ = —ef(x)y and, hence, by substituting it in the last expression, 
an integral equation is obtained for the limit cycle, that is: 



The finite set of limit cycles of Eq. (j5J) also verifies Eq. If the limit cycles are expanded 
in a power series of e, Eq. © imposes the differential relationships that must be verified 




(4) 




(5) 
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for the different orders in e; this is exploited in Sections 3.2 and 3.3. An alternative point 
of view comes from the necessary integral condition ([5]). It can also be exploited to find 
the correct approximate amplitudes of the limit cycles up to a given order in e; this is 
carried out in Section 3.1 up to the order O(e ) and in Section 4 up to any order. 

3 Amplitude and form of the limit cycles 

We proceed now to explain how to calculate recursively the form and amplitude of limit 
cycles of Eq. Q for different orders in e when this parameter is small and f(x) is even. 

3.1 Limit cycles at order zero 

For a given y(x) we define the function S(a) as follows, 

s ( a ) = [ f(x)y(x)dx. (6) 

Then, as it has been established in Eq. © of the previous section, a necessary condition 
for y+(x), hereinafter called y(x), to be a limit cycle of Eq. © with amplitude a is 

S(a) = 0. (7) 

When e — > 0, the limit cycle solutions of Eq. © emerging from the period annulus 
surrounding the center (0, 0) become circles on the plane, 

y (x) = \/ a 2 - x 2 . (8) 

The amplitudes of these limit cycles, a, verify S(a) = 0: 

A)(a) = r f(x)y (x)dx = f f(x)Va 2 -x 2 dx = 0, (9) 

J— a J— a 

and every solution a of Po(a) = is the germ of at least one limit cycle of amplitude a 
|19j . This is a well established result if we recall at this point 0] that /?o( a ) coincides 
with the first Melnikov function or, equivalently, with the Abelian integral defined for the 
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perturbed Hamiltonian system (J2J) whose level curves at e = are given by x 2 + y 2 = a 2 . 
For instance, when f(x) is an even polynomial of degree 2n, (3o( a ) 1S o 2 times a polynomial 
of degree n in a 2 . Then, it has at most n simple positive roots, and we can conclude that 
LMP-conjecture is true in this regime |19j . 



3.2 Limit cycles up to order (9(e ) 



By continuity, we can assume that the limit cycle y(x) for small e is approximated by 
yo(x) plus an expansion in powers of e: 

JV 

y[x) = y (x) + ]T e n y n (x) + 0(e N+1 ). (10) 

n=l 

Replacing this expansion in @ and equating powers of e, every function y n {x) of the 
expansion satisfies a first order linear differential equation of the form: 

(yoVn)' + W n _i = 0, n= 1,2,3,..., N, (11) 

where 

W ,( x ) = l yo( x \f( x ) if n=l, 

The system of equations ()llj) can be iteratively solved for y n {x) by imposing the contour 
condition y n {—o) = V n. Then, we obtain: 



X 



Vn{x) = 1 -~ I W n ^(t)dt 

y (x) J-a 



1 f 1 



yo{x) 



l f x 

-Y,yk{x)y n -k{x)+ y n ^(t)f(t)dt\ . (13) 
1 fc=i J ~ a ) 

Therefore, for small e, the solutions y(x) of Eq. (j2J), which verify y(—a) = 0, are approxi- 
mated up to the order 0{e N ) by 

N 

^-^(ijH^eX^). (14) 

n=0 

For example, up to the order C(e 3 ), t/(x) can be written as 

y (3)( x ) = yo ( x ) + eyi(x) + e 2 y 2 (x) + e 3 y 3 (x), (15) 
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where 



yi(x) = — ^ f y (t)f(t)dt, (16) 

yo{x) J-a 

nix) = -^rs-^f-J X yi(t)f(t)dt, (17) 

2y0{X) yo{X) J-a 

ys(x) = - 4- [ x y2 (t )md t. (18) 

yo{x) yo(x) J-a 



The possible amplitudes a of the limit cycles of Eq. (j3J), up to the order 0(e N ), are 
determined by imposing that the solution y( N \x) vanishes also at x = a: 

N 

y W(a) = 5>"y n (a) = 0. (19) 

n=0 

The condition y( N \a) = identifies the limit cycle y[x) up to the order 0(e N ), but its 
amplitude a only up to the order 0(e N ~ 1 ) because yo(a) vanishes for any value of a. 
Hence, the last condition IJ19JI can be rewritten as 

y W(a)=0 ^e"- 1 y ft (a) = 0. (20) 

n=l 

The solutions a from this last equation are used to calculate y^{x), and hence the approx- 
imation y( N >(x) for the form of the limit cycles can be finally obtained. The iteration of 
this scheme for N = 1, 2, 3 . . . generates the whole expansion of the limit cycle in powers 
of e. 

3.3 Application of the method for the lower orders 
3.3.1 Limit cycles up to the order 0(e 1 ) 



Up to the order 0(1) in the amplitude (N = 1), the above condition (f20|) imposes the 
following restriction over the values of a: 

y (1) (a) = yi(a) = 0. (21) 

The function y\{x) reads 

yi(^) = / yo(t)f(t)dt = _J#,(a)+/ Jto(t)/(t)rftL (22) 

yo{x) J-a yo{x) I Ja J 
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where Po(a) is given in expression @. The function yo(t)f(t) is continuous in the interval 
[—a, a] and of the order 0{y/a — x) when x — * a. Therefore, 

J yo(t)f(t)dt = o({a- x) 3/2 ) when x -> a. (23) 



As yo(x) = Va 2 — x 2 , then, for any a > 0, 

yi ( x ) = -(M®t + 0(a-x), when i ^ a. (24) 
yo{x) 

Hence, from (|22j). a necessary and sufficient condition for to be zero at x = a is just 

the condition Q derived from the integral equation @ for the zero order approximation: 

A) (a) = 0. (25) 

We recover in this way the order O(e ) condition for the amplitudes of the limit cycles 
emerging from the period annulus at order zero ©. The form of these solutions is per- 
turbed at the order 0(e 1 ) and reads: 

y (1 \x) = y (x)+e yi {x), (26) 

with yi (x) given by expression (|16|) when a is solution of Eq. (|25|) . 

3.3.2 Limit cycles up to the order 0(e 2 ) 

Up to order 0(e) in the amplitude (N = 2), condition (|2Uj) imposes the following restriction 
over the values of a: 

y( 2 )(a) = yi(a) + ey 2 (a) = 0. (27) 

Define 

ft (a) = f yi (x)f(x)dx. (28) 

J— a 

Recalling that f(x)yo(x) is an even function of x, and then /?o(a) = 2 J a f(x)yo(x)dx = 
2 f(x)yo(x)dx, we have that: 

ft (a) = - f ^r\dx f f(t)y (t)dt = - f + [ X f(t)y (t)dt) . 

J-ayo{X) J -a J-ayo(X) [2 J ) 
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Because f(t)yo(t) is an even function of t, the t— integral in the above last expression is 
an odd function of x and then vanishes with the exterior integral in x, when it runs in the 
interval [—a, a]. The result is 



A (<0 = -^o(a) f l^dx = -^-(3 (a)(3' (a). 
2 J- a y (x) 2a 

From this equation and Q17|) we find out that 

V2(?) = ^ {\vl(x) - ±/3 (a)f3' (a) + f yx(t)f(t)dt 

yo{x) [2 2a J a 



(29) 



(30) 



Observe that if a is a solution of (3q(o) = 0, that is yi(a) = 0, the same value of a 
automatically satisfies 2/2(0) = 0. This means that, up to the order 0(e), the amplitudes 
of the limit cycles of Eq. (j3J) are the positive roots of the polynomial (3q{o). Then, it seems 
that they do not experience any correction at the order 0(e). Let's show this statement 
more precisely. Write a = ao + a%e + 0(e 2 ) when e — > 0, where ao is a root of Po(a): 
A)( a o) = and 01 is a real number. Then, 



/%(a) = $(a o )aie + 0(e 2 ). 



(31) 



In the foregoing discussion we will consider that the limit e — > must be taken before the 
limit x — > a and then any product of symbols 0(e p ) and ((a — x) q ) can be replaced by 
0(e p ). For any a > 0, yi(x) is a continuous functions in [—a, a] up to the order 0(e) and 
from J21J) and (JSU): 

-A) (a) 



2/1 W 



2/0 (2) 



+ 0(a-x) = 0(e) + 0(a-x). 



(32) 



Then, 



yi (t)f(t)dt = 0((a-x) 2 ) +0(e) 



when x — ► a and e — ► 



(33) 



From the above equation and (|30|) : 

2/2(2;) = j— 

From ((211) and (EU we have that: 



(a-x) 2 +/3i(a) + 0(e) 



(34) 



y( 2 )( X ) = y Q ( X ) + e 



0(a-x) 



/So (a) 



2/0 0) 



+ e z 



0(( 



a — x 



,3/2 



A (a) 



2/o(^) 



+ 0(e 3 



Hence, in order to have y( 2 \a) = up to the order 0(e), the function 

pM(a,e) =/3b(a) + e/3i(o) (35) 

must satisfy 

/3«(a,e) = 0(e 2 ). (36) 

Using that /3 (a) = p' (a ) ai e + 0(e 2 ) and /?i(a) = /?i(a ) + 0(e) = 0(e) (see (fSSjl ). the 
above equation implies 

^(ao)ai = 0. (37) 
Therefore ai = when (3' (ao) ^ 0. The case P' (ao) = 0, where a bifurcation of a multiple 
limit cycle is possible, is not considered here. This means that the amplitude of the limit 
cycles of Eq. © have the form a(e) = ao + a 2 e 2 + C(e 3 ), where a 2 is a real number. This 
is a consequence of the parameter inversion symmetry of Eq. © commented in Section 2. 

The form of these solutions is perturbed and reads: 

y {2) {x) = y (x) + e yi (x) + e 2 y 2 (x), (38) 

where yi(x),y2(x) are given by expressions (|lf)H7j) . with y 2 {x) calculated for the solutions 
a of Po(a) = 0. Observe that as the function yo{x) is an even function of x, then, from 
(|2*2*|) . we see that y±(x) is an odd function of x when a is a root of Po(a), and then, from 
(j!7|) and ()29|) . it is clear that y 2 {x) is an even function of x. 

3.3.3 Limit cycles up to the order (D(e 3 ) 

Up to order 0(e 2 ) in the amplitude (N = 3), condition 1)20(1 imposes the following restric- 
tion over the values of a: 

y(3)(a) = =► yi(a) + ey 2 (a) + e 2 y 3 (a) = 0. (39) 



From the definition (|18| ) of 7/3 (x) we have that 
2/3(«) = 



yo(x) 
1 



yi(x)y 2 (x) + / y 2 (t)f(t)dt 

J — a 

yi{x)y 2 {x) + T y 2 (t)f{t)dt + / y 2 (t)f(t)dt 

J a 
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(40) 
(41) 



From the previous section we have that /3 Q (a) = /3 (a )a 2 e 2 + C(e 3 ) and ft (a) = C(e 2 ). 
Then, from 1)24(1 and ()3fljl we have that, for any a > and up to the order 0(e 2 ), 
and y 2 (x) are continuous functions in [—a, a], yi{x) = 0(a -x) + 0(e 2 ) and 



y 2 (x) = o((a-xfl 2 )+0{e 2 ) 



Then: 



and 



2/3 0) 



1 



— Ua-xfl 2 + ft(a) + 0(e 2 ) , 
y (z) J 

where we have defined 

02(a) = p y 2 {t)f{t)dt. 

J — a 

If the results (|24|). (|^2~)) and (|H|) are substituted in the expression (fTa|) we obtain: 



y (3) W = Vo(x) + e 



0{a-x) 



+e° 



Ua-x)' 



yo(x) 

ft (a 



yo(x) 



+ 0(e 2 



o([a-xf' 2 ) + 0{e 2 ) 
+ 0(e 4 ). 



+ 



Using that /3o(a) = P' (ao)a 2 e 2 + C(e 3 ) and ft(a) = ft(ao) + C( e ) we have that 



y (3) (.x) = 2/o0*0 + e 

0((a-x) 2 ) 



0( a -x)- ^hl e 2 



2/0 0*0 
ft(oo) 



! ((a - x f 2 



+ 



yo{x) 



+ 0(e 4 



(42) 



J X y 2 (t)f{t)dt = O ((a - j;) 5/2 ) + C(e 2 ) when x -» a and e ^ (43) 



(44) 



(45) 



(46) 



(47) 



Hence, in order to have y^\a) = up to the order 0(e 2 ), the function 



/? (2) (a,e) = ft(a) + eft(a) + e 2 ft(a) = [/3 (o )o 2 + ft(o )] e 2 + 0(e 3 ) 



must satisfy 



Therefore 



p( 2 \a,e) = 0(e 3 ). 
ft(oo) 



02 



/?o(«o)' 



The correction of the amplitude up to order 0(e 2 ) is finally obtained: 
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(48) 



(49) 



(50) 



and the form of the solutions at the order 0(e 3 ) reads: 



y 



(x) = y (x) + tyxix) + e 2 y 2 {x) + e 3 y 3 (x) 



(51) 



where yi(x), y2{x), yz(x) are given by expressions I)16I18)> . with y 3 {x) calculated for the 
amplitudes a(e) given by the expression (|50|). 

4 An alternative view: 

the recursion of the integral equation 

Let us see that the method explained in the last section can be viewed in an equivalent 
and alternative way as a recursive approximation in successive powers of e to the integral 
Eq. ©. 

4.1 General algorithm 

Step : Set y (x) = Va 2 - x 2 and N = 0. 

Step 1: Define the approximation of y(x) at order 0{e N ) given by Eq. (|14|): 

N 

y{x)^y( N ){x)^Y.z n yn^)i (52) 

n=0 

and find the solutions a of the integral equation at order 0{e ): 




(53) 



that is equivalent to the equation 



v 



/3( Ar )(a,e)^^6"/3„(a) = 



(54) 



n=0 



where 




(55) 



Step 2: For each solution a of step 1 calculate yN+i(x) by the formula (fT3)l : 



1 




yo(x) 
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i i r 

= rx { ^ E Vk{x)yn-k{x) + / y N (t)f(t)dt . (56) 

Step 3: Replace N by N + 1 and come back to step 1. 

Note: It can be easily found that the application of this algorithm for the lower orders, 
up to N = 2, repeats the results obtained in Section 3.3. 

5 Examples 

The limit cycles in the weak and in the strongly nonlinear regimes of different families 
of Lienard systems were studied in |19j . Here we perform the calculations proposed in 
Section 3 for two concrete examples, which are particular cases of the families 1 and 3 
worked out in 1!) . 

Example 1. The van der Pol oscillator is given for f(x) = x 2 — 1. This system has a 
unique limit cycle, which is stable for e > 0. Hence, the only root of Po(a) is ao = 2. For 
this value of ao, we have y^ 2 '{x) = yo(x) + ey\{x) + t 2 yi{x) with 

y (x) = VA-x 2 , 

and 

y2{x) = 2 -±f{±-x 2 fl 2 . 

We integrate Eq. Q by a Runge-Kutta method in order to obtain the limit cycle. This 
curve is plotted in a continuous trace in Fig. l(a-b) for e = 0.5 and e = 1, respectively. 
The approximated limit cycle y^ 2 \x) is also plotted in those figures with a discontinuous 
trace for the same values of e. Let us remark that, in this case, even up to e = 3, the 
approximation y( 2 \x) to the limit cycle is very good. 

The solution of (a, e) = 0, up to order 0(e 2 ), gives us: 

ft(ao) = Ma) = ^a 2 (a 2 - 4), 
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then 



a(e) = 2 + ±e 2 + 0(e 3 ). (57) 



We must stress at this point that this expansion was also done for the van der Pol system 
in |21| by a perturbation method based on integrating factors. The author reported in 
that article |21| the expansion: 

23 

a(e)=2 + -e 2 + 0(e 3 ). (58) 

Our calculation (|57jl does not agree with this result, but agrees with the computational 
calculation of the 'exact' amplitudes given in Table 1. 

Example 2. The same process is performed for f(x) = 5x 4 — 9x 2 + 1. In this case, 
the system has two limit cycles, one stable and the other unstable. The polynomial (3q{o) 



has two positive roots: ao = y 9 = 0.720677 (unstable limit cycle for e > 0) and 



a = y] = 1.755170 (stable limit cycle for e > 0). 



For the first limit cycle we have y^ 2 \x) = yo(x) + ey\(x) + e 2 y2(x) with 



2/0 (x) = \Jal~x 2 



yi (x) = ^(27 + 3^ - 20x 2 )(x 2 - a 2 ) 



and 



y y ' 2880 



3222 - 218V41 , ^ , , r— , fi 
— + (1003 + 63V41)x 2 - 15(27 + 7V41)x 4 + 200x 6 



The solution of /3^(a, e) = 0, up to order 0(e 2 ), gives us: 

f3 2 (a ) = -0.007357, (5 {a) = — a 2 (5a 4 - 18a 2 + 8). 
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Then 

a(e) = 0.720677 + 0.003908 e 2 + C(e 3 ). (59) 



For the second limit cycle we have y^ 2 \x) = yo(x) + eyi(x) + e 2 y2(x) with 

y (x) = \Jal~ x 2 , 
14 



yi(x) = x 



9 - ViT 5 




and 



V2{x) = 



(ag - x 2 ) 3 / 2 [ 3222 + 218^41 
2880 5 



+ (1003 - 63\/4T)a; 2 + 15(7 VH - 27)x 4 + 200x' 



,6 



The solution of /3^(a, e) = 0, up to order 0(e 2 ), gives us: 



/3 2 (a ) = -0.486199, 



/3 (a) = — a 2 (5a 4 -18a 2 + 8). 



Then 



a(e) = 1.755170 - 0.017880 e 2 + C(e 3 ). 



(60) 



The comparison between the 'exact' and the approximated limit cycles can be seen in 
the plots of Fig. 2(a-b) and Fig. 3(a-b). 

6 Conclusions 

A general algorithm to approximate the form and the amplitude of limit cycles in the 
weakly nonlinear regime of Lienard equations has been presented. In each iteration of the 
method a new order in e is calculated. Thus, the new term added in the expansion of 
e for the form of the limit cycle must verify the differential equation, and the correction 
to this order for the amplitude of the limit cycle is derived from the integral equation. 
Different examples have been worked out and the results are in good agreement with the 
direct integration of the limit cycles. 

In this paper we only have detailed the approximation up to the order e 3 . This process 
can be rewritten in a formalism that allow us to go farther in the order of approximation. 
We are just working in this direction and these new results will be presented elsewhere. 
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Table 1. The value ax represents the approximated amplitude a(e) of the van der Pol 
limit cycle obtained from (|57jl for the indicated values of e. The value cle represents the 
amplitude a obtained by integrating directly the system with a Runge-Kutta method. 



e 


0.1 


0.2 


0.3 


0.4 


0.5 


or 


2.00010 


2.00041 


2.00093 


2.00166 


2.00260 


CLE 


2.00010 


2.00041 


2.00092 


2.00161 


2.00248 


e 


0.6 


0.7 


0.8 


0.9 


1 


ClT 


2.00375 


2.00510 


2.00666 


2.00843 


2.01041 


a E 


2.00351 


2.00466 


2.00591 


2.00724 


2.00862 



Figures 



Fig. la-b: Exact limit cycle (continuous line) and approximated limit cycle (dis- 
continuous line ) up to order C(e 2 ) for the van der Pol system with (a) e = 0.5 and (b) 
e = 1. 

Fig. 2a-b: Exact limit cycle (continuous line) and approximated limit cycle (discon- 
tinuous line ) up to order 0(e 2 ) for the smallest limit cycle of the example 2 with (a) 
e = -0.5 and (b) e = -1. 

Fig. 3a-b: Exact limit cycle (continuous line) and approximated limit cycle (discon- 
tinuous line ) up to order 0{e 2 ) for the biggest limit cycle of the example 2 with (a) e = 0.5 
and (b) e = 1. 
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